package scu.maqiang.homogenization;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES1L21;
import scu.maqiang.fes.FES1L2PE;
import scu.maqiang.mesh.Mesh1L2;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.DirectSSolver;
import scu.maqiang.numeric.DoubleIntegerPair;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.SRMatrix;

public class SOTS1DL2PE {

	public SOTS1DL2PE(FES1L2PE fsPE) {
		this.fsPE = fsPE;
	}

	/**
	 * ͨ������������������ѹ�������ɵ�һά�ṹ�ľ��Ȼ�ѹ�����
	 * 
	 * @param pe1
	 *            ����1��ѹ����ϲ���
	 * @param pe2
	 *            ����2��ѹ����ϲ���
	 * @param ratio
	 *            ����1��ռ�ṹ����, 0��1֮�����
	 */
	public void computeHomogenizedPE(double[] pe1, double[] pe2, double ratio) {
		double G1 = pe1[0];
		double ep1 = pe1[1];
		double e1 = pe1[2];

		double G2 = pe2[0];
		double ep2 = pe2[1];
		double e2 = pe2[2];

		double r1 = ratio, r2 = 1.0 - r1;
		double VG = r1 * G1 / (G1 * ep1 + e1 * e1) + r2 * G2 / (G2 * ep2 + e2 * e2);
		double Vep = r1 * ep1 / (G1 * ep1 + e1 * e1) + r2 * ep2 / (G2 * ep2 + e2 * e2);
		double Ve = r1 * e1 / (G1 * ep1 + e1 * e1) + r2 * e2 / (G2 * ep2 + e2 * e2);

		double V0 = VG * Vep + Ve * Ve;
		System.out.println("G1_0 = " + (VG / V0));
		System.out.println("ep_0 = " + (Vep / V0));
		System.out.println("e1_0 = " + (Ve / V0));
	}

	/**
	 * ͨ������������������ѹ�������ɵ���Գƽṹ�ľ��Ȼ�ѹ�����
	 * 
	 * @param pe1
	 *            ����1��ѹ����ϲ���
	 * @param pe2
	 *            ����2��ѹ����ϲ���
	 * @param ratio
	 *            ����1��ռ�ṹ����, 0��1֮�����
	 */
	public void computeHomogenizedAxisPE(double[] pe1, double[] pe2, double ratio) {
		double G11 = pe1[0];
		double G21 = pe1[1];
		double lambda1 = pe1[2];
		double ep11 = pe1[3];
		double e11 = pe1[4];
		double e21 = pe1[5];

		double G12 = pe2[0];
		double G22 = pe2[1];
		double lambda2 = pe2[2];
		double ep12 = pe2[3];
		double e12 = pe2[4];
		double e22 = pe2[5];

		double r1 = ratio, r2 = 1.0 - r1;
		double VG = r1 * G11 / (G11 * ep11 + e11 * e11) + r2 * G12 / (G12 * ep12 + e12 * e12);
		double Vep = r1 * ep11 / (G11 * ep11 + e11 * e11) + r2 * ep12 / (G12 * ep12 + e12 * e12);
		double Ve = r1 * e11 / (G11 * ep11 + e11 * e11) + r2 * e12 / (G12 * ep12 + e12 * e12);
		double V1 = r1 * (lambda1 * e11 - G11 * e21) / (G11 * ep11 + e11 * e11)
				+ r2 * (lambda2 * e12 - G12 * e22) / (G12 * ep12 + e12 * e12);
		double V2 = r1 * (lambda1 * ep11 + e21 * e11) / (G11 * ep11 + e11 * e11)
				+ r2 * (lambda2 * ep12 + e22 * e12) / (G12 * ep12 + e12 * e12);
		double V0 = VG * Vep + Ve * Ve;
		double G1_0 = VG / V0;
		double ep_0 = Vep / V0;
		double e1_0 = Ve / V0;
		double lambda_0 = e1_0 * V1 + G1_0 * V2;
		double e2_0 = e1_0 * V2 - ep_0 * V1;
		double G2_0 = r1
				* (G21 + ((G11 * e21 - lambda1 * e11) * (e21 - e2_0)
						- (lambda1 * ep11 + e21 * e11) * (lambda1 - lambda_0)) / (G11 * ep11 + e11 * e11))
				+ r2 * (G22 + ((G12 * e22 - lambda2 * e12) * (e22 - e2_0)
						- (lambda2 * ep12 + e22 * e12) * (lambda2 - lambda_0)) / (G12 * ep12 + e12 * e12));
		System.out.println("G1_0 = " + G1_0);
		System.out.println("G2_0 = " + G2_0);
		System.out.println("lambda_0 = " + lambda_0);
		System.out.println("ep_0 = " + ep_0);
		System.out.println("e1_0 = " + e1_0);
		System.out.println("e2_0 = " + e2_0);
	}

	/**
	 * ����һάֱ������ϵ��ѹ�絥�����⼰���Ȼ����ϲ���
	 * 
	 * @param PEFunc
	 *            ѹ��ϵ������, ����ϵ����������, �������6��Ԫ��, ����Ϊѹ�����G1, ep1, e1, �����ܶ�rho, �ṹ�����f1,
	 *            �ṹ����ܶ�rf
	 * @param param
	 *            ѹ�纯����Ӧ����
	 * @param FOCS
	 *            ��ά����, ����Ϊ2, ����Ϊ�����������ɶȸ���, �õ�һ�׵�����
	 * @param D
	 *            2x3 ����, ���صõ��ľ��Ȼ�ϵ���Լ����Ȼ�����������ܶ� D[0][0] = G1_0, D[0][1] = ep1_0,
	 *            D[0][2] = e1_0 D[1][0] = rho_0, D[1][1] = f1_0, D[1][2] = rf_0
	 * @param SOCS
	 *            ��ά����, ����Ϊ3, ����Ϊ�����������ɶȸ������õ����׵�����
	 */
	public void PECellSolution(ScalarFunc[] PEFunc, double[] param, double[][] FOCS, double[][] D,
			double[][] SOCS) {
		SRMatrix A = new SRMatrix(fsPE.GetNdof());
		fsPE.assemblePEStiff(PEFunc, param, BVPType.COMMON, A);
		double[][] RHS = new double[2][fsPE.GetNdof()];
		// ��װһ�׵��������Ҷ���
		fsPE.assembleGlobalVector(PEFunc, param, SOTS1DL2PE::elementPEFirstOrder, BVPType.COMMON, RHS);
		// System.out.println(MVO.toString(RHS[0]));
		// System.out.println(MVO.toString(RHS[1]));
		// ʩ��0�߽�����
		fsPE.applyBC_MBN(A, Direct.All, 1, 2);
		fsPE.applyBC_MBN(RHS[0], Direct.All, 0.0, 1, 2);
		fsPE.applyBC_MBN(RHS[1], Direct.All, 0.0, 1, 2);

		DirectSSolver dsolver = new DirectSSolver(A);
		// NewIterSSolver solver = new NewIterSSolver(A);
		MVO.fill(FOCS, 0.0);
		// ���һ�׵�������
		FOCS[0] = dsolver.LUSolve(RHS[0], 1, 1e-6);
		FOCS[1] = dsolver.LUSolve(RHS[1], 1, 1e-6);
		// solver.CGNR(RHS[0], FOCS[0], 1);
		// solver.CGNR(RHS[1], FOCS[1], 1);
		// System.out.println(MVO.toString(RHS[0]));
		// System.out.println();
		// System.out.println(MVO.toString(RHS[1]));
		// solver.GMRES_SSOR(RHS[0], FOCS[0], 1.50, 20, 1);
		// solver.GMRES_SSOR(RHS[1], FOCS[1], 1.50, 20, 1);

		// ������Ȼ�ϵ��
		computeHomogenizedPE(PEFunc, param, FOCS, D);
		//
		// ��װ���׵��������Ҷ���
		double[][] RHS2 = new double[3][fsPE.GetNdof()];
		fsPE.assembleGlobalVector(PEFunc, param, FOCS, D, SOTS1DL2PE::elementPESecondOrder, BVPType.COMMON, RHS2);

		// ʩ�ӱ߽�����
		MVO.fill(SOCS, 0.0);
		fsPE.applyBC_MBN(RHS2[0], Direct.All, 0.0, 1, 2);
		fsPE.applyBC_MBN(RHS2[1], Direct.All, 0.0, 1, 2);
		fsPE.applyBC_MBN(RHS2[2], Direct.All, 0.0, 1, 2);

		// �����׵�������
		// solver.CGNR(RHS2[0], SOCS[0], 1);
		// solver.CGNR(RHS2[1], SOCS[1], 1);
		// solver.CGNR(RHS2[2], SOCS[2], 1);

		SOCS[0] = dsolver.LUSolve(RHS2[0], 1, 1e-6);
		SOCS[1] = dsolver.LUSolve(RHS2[1], 1, 1e-6);
		SOCS[2] = dsolver.LUSolve(RHS2[2], 1, 1e-6);
	}

	/**
	 * ����һά��Գ�ѹ�絥�����⼰���Ȼ�����ϵ��
	 * 
	 * @param PEFunc
	 *            ѹ��ϵ������, ����ϵ����������, �������9��Ԫ��, ����Ϊѹ�����G1, G2, lambda, ep1, e1, e2
	 *            �����ܶ�rho, �ṹ�����f1, �ṹ����ܶ�rf
	 * @param param
	 *            ѹ�纯����Ӧ����
	 * @param FOCS
	 *            ��ά����, ����Ϊ3, ����Ϊ�����������ɶȸ���, �õ�һ�׵�����
	 * @param D
	 *            3x3 ����, ���صõ��ľ��Ȼ�ϵ���Լ����Ȼ�����������ܶ� D[0][0] = G1_0, D[0][1] = G2_0,
	 *            D[0][2] = lambda_0 D[1][0] = ep1_0, D[1][1] = e1_0, D[1][2] = e2_0
	 *            D[1][0] = rho_0, D[1][1] = f1_0, D[1][2] = rf_0
	 * @param SOCS
	 *            ��ά����, ����Ϊ6, ����Ϊ�����������ɶȸ������õ����׵�����
	 */
	public void AxisPECellSolution(ScalarFunc[] PEFunc, double[] param, double[][] FOCS, double[][] D,
			double[][] SOCS) {
		SRMatrix A = new SRMatrix(fsPE.GetNdof());
		ScalarFunc[] cellPEFunc = new ScalarFunc[] { PEFunc[0], PEFunc[3], PEFunc[4] };
		fsPE.assemblePEStiff(cellPEFunc, param, BVPType.COMMON, A);
		double[][] RHS = new double[3][fsPE.GetNdof()];
		fsPE.assembleGlobalVector(PEFunc, param, SOTS1DL2PE::elementAxisPEFirstOrder, BVPType.COMMON, RHS);
		// System.out.println(MVO.toString(RHS[0]));
		// System.out.println(MVO.toString(RHS[1]));
		// System.out.println(MVO.toString(RHS[2]));
		fsPE.applyBC_MBN(A, Direct.All, 1, 2);
		fsPE.applyBC_MBN(RHS[0], Direct.All, 0.0, 1, 2);
		fsPE.applyBC_MBN(RHS[1], Direct.All, 0.0, 1, 2);
		fsPE.applyBC_MBN(RHS[2], Direct.All, 0.0, 1, 2);

		DirectSSolver dsolver = new DirectSSolver(A);
		// NewIterSSolver solver = new NewIterSSolver(A);
		MVO.fill(FOCS, 0.0);
		FOCS[0] = dsolver.LUSolve(RHS[0], 1, 1e-6);
		FOCS[1] = dsolver.LUSolve(RHS[1], 1, 1e-6);
		FOCS[2] = dsolver.LUSolve(RHS[2], 1, 1e-6);
		//
		computeHomogenizedAxisPE(PEFunc, param, FOCS, D);
		//
		double[][] RHS2 = new double[6][fsPE.GetNdof()];
		fsPE.assembleGlobalVector(PEFunc, param, FOCS, D, SOTS1DL2PE::elementAxisPESecondOrder, BVPType.COMMON,
				RHS2);

		MVO.fill(SOCS, 0.0);
		fsPE.applyBC_MBN(RHS2[0], Direct.All, 0.0, 1, 2);
		fsPE.applyBC_MBN(RHS2[1], Direct.All, 0.0, 1, 2);
		fsPE.applyBC_MBN(RHS2[2], Direct.All, 0.0, 1, 2);
		fsPE.applyBC_MBN(RHS2[3], Direct.All, 0.0, 1, 2);
		fsPE.applyBC_MBN(RHS2[4], Direct.All, 0.0, 1, 2);
		fsPE.applyBC_MBN(RHS2[5], Direct.All, 0.0, 1, 2);
		for(int i = 0; i < 6; i++) {
			SOCS[i] = dsolver.LUSolve(RHS2[0], 0, 1e-6);
		}
	}

	public void SphericalPECellSolution(ScalarFunc[] PEFunc, double[] param, double[][] FOCS, double[][] D,
			double[][] SOCS) {

	}

	/**
	 * ѹ��һ�׵������ⵥԪ�Ҷ�����
	 * 
	 * @param xy
	 *            ��Ԫ�������
	 * @param coef
	 *            ��Ԫ����ϵ��
	 * @param tp
	 *            ��������
	 * @param eleXY
	 *            ��Ԫ�Ҷ�����
	 */
	public static void elementPEFirstOrder(int[] ele, double[][] xy, double[][] coef, BVPType tp, double[][] eleXY) {
		double mG = 0.5 * (coef[0][0] + coef[0][1]);
		double mEp = 0.5 * (coef[1][0] + coef[1][1]);
		double mE1 = 0.5 * (coef[2][0] + coef[2][1]);
		eleXY[0][0] = mG;
		eleXY[0][1] = -mE1;
		eleXY[0][2] = -mG;
		eleXY[0][3] = mE1;
		eleXY[1][0] = mE1;
		eleXY[1][1] = mEp;
		eleXY[1][2] = -mE1;
		eleXY[1][3] = -mEp;
	}

	/**
	 * һά��Գ�ѹ�絥�������Ҷ�����
	 * 
	 * @param xy
	 *            ��Ԫ����
	 * @param coef
	 *            ��Գ�ѹ����ϲ���
	 * @param tp
	 *            ��������
	 * @param eleXY
	 *            ��Ԫ�Ҷ�����
	 */
	public static void elementAxisPEFirstOrder(int[] ele, double[][] xy, double[][] coef, BVPType tp, double[][] eleXY) {
		double mG1 = 0.5 * (coef[0][0] + coef[0][1]);
		// double mG2 = 0.5 * (coef[1][0] + coef[1][1]);
		double lambda = 0.5 * (coef[2][0] + coef[2][1]);
		double mEp = 0.5 * (coef[3][0] + coef[3][1]);
		double mE1 = 0.5 * (coef[4][0] + coef[4][1]);
		double mE2 = 0.5 * (coef[5][1] + coef[5][1]);
		eleXY[0][0] = mG1;
		eleXY[0][1] = -mE1;
		eleXY[0][2] = -mG1;
		eleXY[0][3] = mE1;

		eleXY[1][0] = mE1;
		eleXY[1][1] = mEp;
		eleXY[1][2] = -mE1;
		eleXY[1][3] = -mEp;

		eleXY[2][0] = lambda;
		eleXY[2][1] = -mE2;
		eleXY[2][2] = -lambda;
		eleXY[2][3] = mE2;
	}

	/**
	 * 
	 * @param kxxFunc
	 * @param param
	 * @param FOCS
	 * @param D
	 *            2x3 ����, ���صõ��ľ��Ȼ�ϵ���Լ����Ȼ�����������ܶ� D[0][0] = G1_0, D[0][1] = ep1_0,
	 *            D[0][2] = e1_0 D[1][0] = rho_0, D[1][1] = f1_0, D[1][2] = rf_0
	 */
	public void computeHomogenizedPE(ScalarFunc[] kxxFunc, double[] param, double[][] FOCS, double[][] D) {
		double[] dx = new double[2];
		double[][] xy = new double[2][1];
		Mesh1L2 mesh = (Mesh1L2) fsPE.getMesh();
		double[] N1H = new double[4];
		double[] FR1 = new double[4];
		double N1Dx, HDx, R1Dx, FDx;
		double allLength = 0.0;
		double[][] coef = new double[kxxFunc.length][2];
		double length = 0.0;
		MVO.fill(D, 0.0);
		for (int i = 0; i < mesh.getNt(); i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, xy);
			int[] dof = fsPE.getDof(ele);
			MVO.getSlice(FOCS[0], dof, N1H);
			MVO.getSlice(FOCS[1], dof, FR1);
			double[] N1 = new double[] { N1H[0], N1H[2] };
			double[] H = new double[] { N1H[1], N1H[3] };
			double[] F = new double[] { FR1[0], FR1[2] };
			double[] R1 = new double[] { FR1[1], FR1[3] };
			length = Math.abs(xy[0][0] - xy[1][0]);
			dx[0] = -1.0 / length;
			dx[1] = 1.0 / length;

			N1Dx = MVO.dot_product(N1, dx);
			HDx = MVO.dot_product(H, dx);
			R1Dx = MVO.dot_product(R1, dx);
			FDx = MVO.dot_product(F, dx);
			allLength += length;
			fsPE.getCoefFromFunc(kxxFunc, xy, mesh.getElementLabel(i), param, coef);
			double mG = 0.5 * (coef[0][0] + coef[0][1]);
			double mEp = 0.5 * (coef[1][0] + coef[1][1]);
			double mE1 = 0.5 * (coef[2][0] + coef[2][1]);
			double mRho = 0.5 * (coef[3][0] + coef[3][1]);
			double mf = 0.5 * (coef[4][0] + coef[4][1]);
			double mrf = 0.5 * (coef[5][0] + coef[5][1]);

			D[0][0] += (mG * (1 + N1Dx) + mE1 * HDx) * length; // G1_0
			D[0][1] += (mEp * (1 + R1Dx) - mE1 * FDx) * length; // ep1_0
			D[0][2] += (mE1 * (1 + R1Dx) + mG * FDx) * length; // e1_0
			D[1][0] += mRho * length; // rho_0
			D[1][1] += mf * length; // f1_0
			D[1][2] += mrf * length; // rf_0
		}
		MVO.mulSelf(D, allLength);
	}

	/**
	 * ������Գ�ѹ����Ͼ��Ȼ�����
	 * 
	 * @param kxxFunc
	 * @param param
	 * @param FOCS
	 * @param D
	 */
	public void computeHomogenizedAxisPE(ScalarFunc[] kxxFunc, double[] param, double[][] FOCS, double[][] D) {
		double[] dx = new double[2];
		double[][] xy = new double[2][1];
		Mesh1L2 mesh = (Mesh1L2) fsPE.getMesh();
		double[] N1H = new double[4];
		double[] FR1 = new double[4];
		double[] MS = new double[4];
		double N1Dx, HDx, R1Dx, FDx, MDx, SDx;
		double allLength = 0.0;
		double[][] coef = new double[kxxFunc.length][2];
		double length = 0.0;
		MVO.fill(D, 0.0);
		for (int i = 0; i < mesh.getNt(); i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, xy);
			int[] dof = fsPE.getDof(ele);
			MVO.getSlice(FOCS[0], dof, N1H);
			MVO.getSlice(FOCS[1], dof, FR1);
			MVO.getSlice(FOCS[2], dof, MS);
			double[] N1 = new double[] { N1H[0], N1H[2] };
			double[] H = new double[] { N1H[1], N1H[3] };
			double[] F = new double[] { FR1[0], FR1[2] };
			double[] R1 = new double[] { FR1[1], FR1[3] };
			double[] M = new double[] { MS[0], MS[2] };
			double[] S = new double[] { MS[1], MS[3] };
			length = Math.abs(xy[0][0] - xy[1][0]);
			dx[0] = -1.0 / length;
			dx[1] = 1.0 / length;

			N1Dx = MVO.dot_product(N1, dx);
			HDx = MVO.dot_product(H, dx);
			R1Dx = MVO.dot_product(R1, dx);
			FDx = MVO.dot_product(F, dx);
			MDx = MVO.dot_product(M, dx);
			SDx = MVO.dot_product(S, dx);
			allLength += length;
			fsPE.getCoefFromFunc(kxxFunc, xy, mesh.getElementLabel(i), param, coef);
			double mG1 = 0.5 * (coef[0][0] + coef[0][1]);
			double mG2 = 0.5 * (coef[1][0] + coef[1][1]);
			double mLambda = 0.5 * (coef[2][0] + coef[2][1]);
			double mEp = 0.5 * (coef[3][0] + coef[3][1]);
			double mE1 = 0.5 * (coef[4][0] + coef[4][1]);
			double mE2 = 0.5 * (coef[5][1] + coef[5][1]);
			double mRho = 0.5 * (coef[6][0] + coef[6][1]);
			double mf = 0.5 * (coef[7][0] + coef[7][1]);
			double mrf = 0.5 * (coef[8][0] + coef[8][1]);

			D[0][0] += (mG1 * (1 + N1Dx) + mE1 * HDx) * length; // G1_0
			D[0][1] += (mG2 + mLambda * MDx + mE2 * SDx) * length; // G2_0
			D[0][2] += (mLambda * (1 + N1Dx) + mE2 * HDx) * length; // lambda_0
			D[1][0] += (mEp * (1 + R1Dx) - mE1 * FDx) * length; // ep1_0
			D[1][1] += (mE1 * (1 + R1Dx) + mG1 * FDx) * length; // e1_0
			D[1][2] += (mE2 * (1 + R1Dx) + mLambda * FDx) * length; // e2_0
			D[2][0] += mRho * length; // rho_0
			D[2][1] += mf * length; // f1_0
			D[2][2] += mrf * length; // rf_0
		}
		MVO.mulSelf(D, allLength);
	}

	/**
	 * 
	 * @param coord
	 * @param coef
	 * @param FOCS
	 * @param D
	 * @param tp
	 * @param eleXY
	 */
	public static void elementPESecondOrder(int[] ele, double[][] coord, double[][] coef, double[][] FOCS, double[][] D,
			BVPType tp, double[][] eleXY) {
		double[] x = MVO.col(coord, 0);
		double L = Math.abs(x[1] - x[0]);
		double[] dx = { -1.0 / L, 1.0 / L };
		double[] N1H = FOCS[0];
		double[] FR1 = FOCS[1];
		// double[] N1 = new double[] {N1H[0], N1H[2]};
		// double[] H = new double[] {N1H[1], N1H[3]};
		// double[] R1 = new double[] {R1F[0], R1F[2]};
		// double[] F = new double[] {R1F[1], R1F[3]};
		//
		// N1Dx = MVO.dot_product(N1, dx);
		// HDx = MVO.dot_product(H, dx);
		// R1Dx = MVO.dot_product(R1, dx);
		// FDx = MVO.dot_product(F, dx);

		double G1_0 = D[0][0];
		double ep1_0 = D[0][1];
		double e1_0 = D[0][2];
		double rho_0 = D[1][0];
		double f1_0 = D[1][1];
		double rf_0 = D[1][2];

		double[] G1 = coef[0];
		double[] ep1 = coef[1];
		double[] e1 = coef[2];
		double[] rho = coef[3];
		double[] f1 = coef[4];
		double[] rf = coef[5];
		// double mG = 0.5 * (coef[0][0] + coef[0][1]);
		// double mEp = 0.5 * (coef[1][0] + coef[1][1]);
		// double mE1 = 0.5 * (coef[2][0] + coef[2][1]);
		// double mRho = 0.5 * (coef[3][0] + coef[3][1]);
		// double mf = 0.5 * (coef[4][0] + coef[4][1]);
		// double mrf = 0.5 * (coef[5][0] + coef[5][1]);

		double mG1_1 = (2 * G1[0] + G1[1]) / 6.0;
		double mG1_2 = (G1[0] + 2 * G1[1]) / 6.0;
		double mEp1_1 = (2 * ep1[0] + ep1[1]) / 6.0;
		double mEp1_2 = (ep1[0] + 2 * ep1[1]) / 6.0;
		double mE1_1 = (2 * e1[0] + e1[1]) / 6.0;
		double mE1_2 = (e1[0] + 2 * e1[1]) / 6.0;
		double mRho_1 = (2 * rho[0] + rho[1]) / 6.0;
		double mRho_2 = (rho[0] + 2 * rho[1]) / 6.0;
		// double mf1_1 = (2 * f1[0] + f1[1]) / 6.0;
		// double mf1_2 = (f1[0] + 2 * f1[1]) / 6.0;
		double mrf_1 = (2 * rf[0] + rf[1]) / 6.0;
		double mrf_2 = (rf[0] + 2 * rf[1]) / 6.0;

		double[] Ge = new double[] { mG1_1, mE1_1, mG1_2, mE1_2 };
		double[] eep = new double[] { -mE1_1, mEp1_1, -mE1_2, mEp1_2 };
		double dN1H_1 = MVO.dot_product(Ge, N1H);
		double dN1H_2 = MVO.dot_product(eep, N1H);
		eleXY[0][0] = dN1H_1 + G1_0 * (0.5 - mRho_1 / rho_0) * L;
		eleXY[0][1] = dN1H_2 - e1_0 * (0.5 - mrf_1 / rf_0) * L;
		eleXY[0][2] = -dN1H_1 + G1_0 * (0.5 - mRho_2 / rho_0) * L;
		eleXY[0][3] = -dN1H_2 - e1_0 * (0.5 - mrf_2 / rf_0) * L;

		double dFR1_1 = MVO.dot_product(Ge, FR1);
		double dFR1_2 = MVO.dot_product(eep, FR1);
		eleXY[1][0] = dFR1_1 + e1_0 * (0.5 - mRho_1 / rho_0) * L;
		eleXY[1][1] = dFR1_2 + ep1_0 * (0.5 - mrf_1 / rf_0) * L;
		eleXY[1][2] = -dFR1_1 + e1_0 * (0.5 - mRho_2 / rho_0) * L;
		eleXY[1][3] = -dFR1_2 + ep1_0 * (0.5 - mrf_2 / rf_0) * L;

		eleXY[2][0] = (mrf_1 - f1_0 * mRho_1 / rho_0) * L;
		eleXY[2][2] = (mrf_2 - f1_0 * mRho_2 / rho_0) * L;
	}

	/**
	 * 
	 * @param coord
	 * @param coef
	 * @param FOCS
	 * @param D
	 * @param tp
	 * @param eleXY
	 */
	public static void elementAxisPESecondOrder(int[] ele, double[][] coord, double[][] coef, double[][] FOCS, double[][] D,
			BVPType tp, double[][] eleXY) {
		double[] x = MVO.col(coord, 0);
		double L = Math.abs(x[1] - x[0]);
		double[] dx = { -1.0 / L, 1.0 / L };
		double[] N1H = FOCS[0];
		double[] FR1 = FOCS[1];
		double[] MS = FOCS[2];
		// double[] N1 = new double[] {N1H[0], N1H[2]};
		// double[] H = new double[] {N1H[1], N1H[3]};
		// double[] R1 = new double[] {R1F[0], R1F[2]};
		// double[] F = new double[] {R1F[1], R1F[3]};
		//
		// N1Dx = MVO.dot_product(N1, dx);
		// HDx = MVO.dot_product(H, dx);
		// R1Dx = MVO.dot_product(R1, dx);
		// FDx = MVO.dot_product(F, dx);

		double G1_0 = D[0][0];
		double G2_0 = D[0][1];
		double lambda_0 = D[0][2];
		double ep1_0 = D[1][0];
		double e1_0 = D[1][1];
		double e2_0 = D[1][2];
		double rho_0 = D[2][0];
		double f1_0 = D[2][1];
		double rf_0 = D[2][2];

		double[] G1 = coef[0];
		double[] G2 = coef[1];
		double[] lambda = coef[2];
		double[] ep1 = coef[3];
		double[] e1 = coef[4];
		double[] e2 = coef[5];
		double[] rho = coef[6];
		double[] f1 = coef[7];
		double[] rf = coef[8];
		// double mG = 0.5 * (coef[0][0] + coef[0][1]);
		// double mEp = 0.5 * (coef[1][0] + coef[1][1]);
		// double mE1 = 0.5 * (coef[2][0] + coef[2][1]);
		// double mRho = 0.5 * (coef[3][0] + coef[3][1]);
		// double mf = 0.5 * (coef[4][0] + coef[4][1]);
		// double mrf = 0.5 * (coef[5][0] + coef[5][1]);

		double mG1_1 = (2 * G1[0] + G1[1]) / 6.0;
		double mG1_2 = (G1[0] + 2 * G1[1]) / 6.0;
		double mEp1_1 = (2 * ep1[0] + ep1[1]) / 6.0;
		double mEp1_2 = (ep1[0] + 2 * ep1[1]) / 6.0;
		double mE1_1 = (2 * e1[0] + e1[1]) / 6.0;
		double mE1_2 = (e1[0] + 2 * e1[1]) / 6.0;
		double mRho_1 = (2 * rho[0] + rho[1]) / 6.0;
		double mRho_2 = (rho[0] + 2 * rho[1]) / 6.0;
		// double mf1_1 = (2 * f1[0] + f1[1]) / 6.0;
		// double mf1_2 = (f1[0] + 2 * f1[1]) / 6.0;
		double mrf_1 = (2 * rf[0] + rf[1]) / 6.0;
		double mrf_2 = (rf[0] + 2 * rf[1]) / 6.0;

		// N11H11
		double[] Ge = new double[] { mG1_1, mE1_1, mG1_2, mE1_2 };
		double[] eep = new double[] { -mE1_1, mEp1_1, -mE1_2, mEp1_2 };
		double dN1H_1 = MVO.dot_product(Ge, N1H);
		double dN1H_2 = MVO.dot_product(eep, N1H);
		eleXY[0][0] = dN1H_1 + G1_0 * (0.5 - mRho_1 / rho_0) * L;
		eleXY[0][1] = dN1H_2 + e1_0 * (0.5 - mrf_1 / rf_0) * L;
		eleXY[0][2] = -dN1H_1 + G1_0 * (0.5 - mRho_2 / rho_0) * L;
		eleXY[0][3] = -dN1H_2 + e1_0 * (0.5 - mrf_2 / rf_0) * L;

		// F11R11
		double dFR1_1 = MVO.dot_product(Ge, FR1);
		double dFR1_2 = MVO.dot_product(eep, FR1);
		eleXY[1][0] = dFR1_1 + e1_0 * (0.5 - mRho_1 / rho_0) * L;
		eleXY[1][1] = dFR1_2 - ep1_0 * (0.5 - mrf_1 / rf_0) * L;
		eleXY[1][2] = -dFR1_1 + e1_0 * (0.5 - mRho_2 / rho_0) * L;
		eleXY[1][3] = -dFR1_2 - ep1_0 * (0.5 - mrf_2 / rf_0) * L;

		// PW
		eleXY[2][0] = (mrf_1 - f1_0 * mRho_1 / rho_0) * L;
		eleXY[2][2] = (mrf_2 - f1_0 * mRho_2 / rho_0) * L;

		////
		// NN1HH1
		double dMS_1 = MVO.dot_product(Ge, MS);
		double dMS_2 = MVO.dot_product(eep, MS);
		double[] N1 = new double[] { N1H[0], N1H[2] };
		double[] H = new double[] { N1H[1], N1H[3] };
		double N1lambda = (MVO.sum(lambda) * MVO.sum(N1) + MVO.dot_product(N1, lambda)) / 6.0;
		double N1E2 = (MVO.sum(e2) * MVO.sum(N1) + MVO.dot_product(e2, N1)) / 6.0;

		double mLambda = MVO.average(lambda);
		double mE2 = MVO.average(e2);
		double mRho = MVO.average(rho);
		double mRf = MVO.average(rf);
		double temp1 = (mLambda * MVO.dot_product(N1, dx) + mE2 * MVO.dot_product(H, dx) + mLambda - lambda_0) * 0.5
				* L;
		double temp2 = 0.5 * (e1_0 + e2_0) * (1 - mRf / rf_0) * L;
		eleXY[3][0] = dMS_1 + N1lambda - temp1;
		eleXY[3][1] = dMS_2 - N1E2 - temp2;
		eleXY[3][2] = -dMS_1 - N1lambda - temp1;
		eleXY[3][3] = -dMS_2 + N1E2 - temp2;

		// MM1SS1
		double[] M = new double[] { MS[0], MS[2] };
		double[] S = new double[] { MS[1], MS[3] };
		double mG2 = MVO.average(G2);
		double Mlambda = (MVO.sum(lambda) * MVO.sum(M) + MVO.dot_product(M, lambda)) / 6.0;
		double ME2 = (MVO.sum(e2) * MVO.sum(M) + MVO.dot_product(e2, M)) / 6.0;
		temp1 = (mLambda * MVO.dot_product(M, dx) + mE2 * MVO.dot_product(S, dx) + mG2 - mRho / rho_0 * G2_0) * 0.5 * L;
		eleXY[4][0] = -dMS_1 + Mlambda - temp1;
		eleXY[4][1] = -dMS_2 - ME2;
		eleXY[4][2] = dMS_1 - Mlambda - temp1;
		eleXY[4][3] = dMS_2 + ME2;

		// FF1RR1
		double[] F = new double[] { FR1[0], FR1[2] };
		double[] R1 = new double[] { FR1[1], FR1[3] };
		double mEp = MVO.average(ep1);

		double Flambda = (MVO.sum(lambda) * MVO.sum(F) + MVO.dot_product(F, lambda)) / 6.0;
		double FE2 = (MVO.sum(e2) * MVO.sum(F) + MVO.dot_product(e2, F)) / 6.0;
		temp1 = (mLambda * MVO.dot_product(F, dx) + mE2 * MVO.dot_product(R1, dx) + mE2 - mRho * e2_0 / rho_0
				- e1_0 * (1 - mRho / rho_0)) * 0.5 * L;
		temp2 = (mEp * (1 - mRf / rf_0)) * 0.5 * L;
		eleXY[5][0] = Flambda - temp1;
		eleXY[5][1] = -FE2 + temp2;
		eleXY[5][2] = -Flambda - temp1;
		eleXY[5][3] = FE2 + temp2;
	}

	/**
	 * ��װһάѹ����Ͼ��Ȼ���һ������׽��ƽ�
	 * 
	 * @param fsFine
	 * @param T_0
	 * @param T_1
	 * @param T_2
	 * @param fsHomo
	 * @param T_Homo
	 * @param fsCell
	 * @param FOCS
	 * @param SOCS
	 * @param N
	 * @param ep
	 */
	public void assemblePESOTSSolution(FES1L2PE fsFine, double[] T_0, double[] T_1, double[] T_2, FES1L2PE fsHomo,
                                       double[] T_Homo, FES1L2PE fsCell, double[][] FOCS, double[][] SOCS, int N, double ep) {
		double[] xy = new double[1];
		double[] cellXY = new double[1];

		double[] eleValHomo = new double[2];
		double[] eleValCell = new double[2];

		FES1L21 fsHH = new FES1L21((Mesh1L2) (fsHomo.mesh));

		int homoNv = fsHomo.mesh.getNv();
		double[][] UPhi = new double[2][homoNv];
		fsHomo.extractComponent(T_Homo, UPhi);

		double[][] gradU = new double[1][homoNv];
		double[][] gradPhi = new double[1][homoNv];
		double[][] SeU = new double[1][homoNv];
		double[][] SePhi = new double[1][homoNv];

		fsHH.computeGradient(UPhi[0], gradU);
		fsHH.computeGradient(UPhi[1], gradPhi);
		fsHH.computeGradient(gradU[0], SeU);
		fsHH.computeGradient(gradPhi[0], SePhi);

		FES1L21 fsCC = new FES1L21((Mesh1L2) (fsCell.mesh));
		int cellNv = fsCell.mesh.getNv();
		double[][] N1H = new double[2][cellNv];
		double[][] FR1 = new double[2][cellNv];
		fsCell.extractComponent(FOCS[0], N1H);
		fsCell.extractComponent(FOCS[1], FR1);

		double[][] N11H11 = new double[2][cellNv];
		double[][] F11R11 = new double[2][cellNv];
		double[][] PW = new double[2][cellNv];
		fsCell.extractComponent(SOCS[0], N11H11);
		fsCell.extractComponent(SOCS[1], F11R11);
		fsCell.extractComponent(SOCS[2], PW);

		double corrector11, corrector12, corrector21, corrector22;

		DoubleIntegerPair homoData = new DoubleIntegerPair(0.0, 0);
		int[] eleIdxHomo = fsHomo.mesh.getElement(0).clone();

		DoubleIntegerPair cellData = new DoubleIntegerPair(0.0, 0);
		int[] eleIdxCell = fsCell.mesh.getElement(0).clone();
		double val = 0.0;
		for (int i = 0, n = fsFine.mesh.getNv(); i < n; i++) {
			fsFine.mesh.getCoordPerNode(i, xy);
			cellXY[0] = xy[0] * N - (int) (xy[0] * N);
			fsHH.getValue(UPhi[0], xy, homoData, eleIdxHomo, eleValHomo);
			T_0[2 * i] = homoData.first;
			T_0[2 * i + 1] = MVO.dot_product(UPhi[1], eleIdxHomo, eleValHomo);

			fsCC.getValue(N1H[0], cellXY, cellData, eleIdxCell, eleValCell);
			val = cellData.first;

			corrector11 = val * MVO.dot_product(gradU[0], eleIdxHomo, eleValHomo);
			corrector11 += MVO.dot_product(FR1[0], eleIdxCell, eleValCell)
					* MVO.dot_product(gradPhi[0], eleIdxHomo, eleValHomo);
			T_1[2 * i] = T_0[2 * i] + ep * corrector11;
			corrector12 = MVO.dot_product(FR1[1], eleIdxCell, eleValCell)
					* MVO.dot_product(gradPhi[0], eleIdxHomo, eleValHomo);
			corrector12 += MVO.dot_product(N1H[1], eleIdxCell, eleValCell)
					* MVO.dot_product(gradU[0], eleIdxHomo, eleValHomo);
			T_1[2 * i + 1] = T_0[2 * i + 1] + ep * corrector12;

			corrector21 = MVO.dot_product(N11H11[0], eleIdxCell, eleValCell)
					* MVO.dot_product(SeU[0], eleIdxHomo, eleValHomo);
			corrector21 += MVO.dot_product(F11R11[0], eleIdxCell, eleValCell)
					* MVO.dot_product(SePhi[0], eleIdxHomo, eleValHomo);
			corrector21 += MVO.dot_product(PW[0], eleIdxCell, eleValCell);
			T_2[2 * i] = T_1[2 * i] + ep * ep * corrector21;

			corrector22 = MVO.dot_product(N11H11[1], eleIdxCell, eleValCell)
					* MVO.dot_product(SeU[0], eleIdxHomo, eleValHomo);
			corrector22 += MVO.dot_product(F11R11[1], eleIdxCell, eleValCell)
					* MVO.dot_product(SePhi[0], eleIdxHomo, eleValHomo);
			corrector22 += MVO.dot_product(PW[1], eleIdxCell, eleValCell);
			T_2[2 * i + 1] = T_1[2 * i + 1] + ep * ep * corrector22;
		}
	}

	/**
	 * ��װƽ����Գ�ѹ����Ͼ��Ȼ���һ������׽��ƽ�
	 * 
	 * @param fsFine
	 * @param T_0
	 * @param T_1
	 * @param T_2
	 * @param fsHomo
	 * @param T_Homo
	 * @param fsCell
	 * @param FOCS
	 * @param SOCS
	 * @param N
	 * @param ep
	 */
	public void assembleAxisPESOTSSolution(FES1L2PE fsFine, double[] T_0, double[] T_1, double[] T_2, FES1L2PE fsHomo,
                                           double[] T_Homo, FES1L2PE fsCell, double[][] FOCS, double[][] SOCS, int N, double ep) {
		double[] xy = new double[1];
		double[] cellXY = new double[1];

		double[] eleValHomo = new double[2];
		double[] eleValCell = new double[2];

		FES1L21 fsHH = new FES1L21((Mesh1L2) (fsHomo.mesh));

		int homoNv = fsHomo.mesh.getNv();
		double[][] UPhi = new double[2][homoNv];
		fsHomo.extractComponent(T_Homo, UPhi);

		double[][] gradU = new double[1][homoNv];
		double[][] gradPhi = new double[1][homoNv];
		double[][] SeU = new double[1][homoNv];
		double[][] SePhi = new double[1][homoNv];

		fsHH.computeGradient(UPhi[0], gradU);
		fsHH.computeGradient(UPhi[1], gradPhi);
		fsHH.computeGradient(gradU[0], SeU);
		fsHH.computeGradient(gradPhi[0], SePhi);

		FES1L21 fsCC = new FES1L21((Mesh1L2) (fsCell.mesh));
		int cellNv = fsCell.mesh.getNv();
		double[][] N1H = new double[2][cellNv];
		double[][] FR1 = new double[2][cellNv];
		double[][] MS = new double[2][cellNv];
		fsCell.extractComponent(FOCS[0], N1H);
		fsCell.extractComponent(FOCS[1], FR1);
		fsCell.extractComponent(FOCS[2], MS);

		double[][] N11H11 = new double[2][cellNv];
		double[][] F11R11 = new double[2][cellNv];
		double[][] PW = new double[2][cellNv];
		double[][] NN1HH1 = new double[2][cellNv];
		double[][] MM1SS1 = new double[2][cellNv];
		double[][] FF1RR1 = new double[2][cellNv];
		fsCell.extractComponent(SOCS[0], N11H11);
		fsCell.extractComponent(SOCS[1], F11R11);
		fsCell.extractComponent(SOCS[2], PW);
		fsCell.extractComponent(SOCS[3], NN1HH1);
		fsCell.extractComponent(SOCS[4], MM1SS1);
		fsCell.extractComponent(SOCS[5], FF1RR1);

		double corrector11, corrector12, corrector21, corrector22;

		DoubleIntegerPair homoData = new DoubleIntegerPair(0.0, 0);
		int[] eleIdxHomo = fsHomo.mesh.getElement(0).clone();

		DoubleIntegerPair cellData = new DoubleIntegerPair(0.0, 0);
		int[] eleIdxCell = fsCell.mesh.getElement(0).clone();
		double val = 0.0;
		double uVal = 0.0;
		for (int i = 0, n = fsFine.mesh.getNv(); i < n; i++) {
			fsFine.mesh.getCoordPerNode(i, xy);
			cellXY[0] = xy[0] * N - (int) (xy[0] * N);
			//������������⣡������������������������������������
			fsHH.getValue(UPhi[0], xy, homoData, eleIdxHomo, eleValHomo);
			uVal = homoData.first;
			T_0[2 * i] = homoData.first;
			T_0[2 * i + 1] = MVO.dot_product(UPhi[1], eleIdxHomo, eleValHomo);

			fsCC.getValue(N1H[0], cellXY, cellData, eleIdxCell, eleValCell);
			val = cellData.first;

			corrector11 = val * MVO.dot_product(gradU[0], eleIdxHomo, eleValHomo);
			corrector11 += MVO.dot_product(FR1[0], eleIdxCell, eleValCell)
					* MVO.dot_product(gradPhi[0], eleIdxHomo, eleValHomo);
			corrector11 += MVO.dot_product(MS[0], eleIdxCell, eleValCell) * uVal / xy[0];
			T_1[2 * i] = T_0[2 * i] + ep * corrector11;

			corrector12 = MVO.dot_product(FR1[1], eleIdxCell, eleValCell)
					* MVO.dot_product(gradPhi[0], eleIdxHomo, eleValHomo);
			corrector12 += MVO.dot_product(N1H[1], eleIdxCell, eleValCell)
					* MVO.dot_product(gradU[0], eleIdxHomo, eleValHomo);
			corrector12 += MVO.dot_product(MS[1], eleIdxCell, eleValCell) * uVal / xy[0];
			T_1[2 * i + 1] = T_0[2 * i + 1] + ep * corrector12;

			corrector21 = MVO.dot_product(N11H11[0], eleIdxCell, eleValCell)
					* MVO.dot_product(SeU[0], eleIdxHomo, eleValHomo);
			corrector21 += MVO.dot_product(F11R11[0], eleIdxCell, eleValCell)
					* MVO.dot_product(SePhi[0], eleIdxHomo, eleValHomo);
			corrector21 += MVO.dot_product(PW[0], eleIdxCell, eleValCell);
			corrector21 += MVO.dot_product(NN1HH1[0], eleIdxCell, eleValCell)
					* MVO.dot_product(gradU[0], eleIdxHomo, eleValHomo) / xy[0];
			corrector21 += MVO.dot_product(MM1SS1[0], eleIdxCell, eleValCell) * uVal / (xy[0] * xy[0]);
			corrector21 += MVO.dot_product(FF1RR1[0], eleIdxCell, eleValCell)
					* MVO.dot_product(gradPhi[0], eleIdxHomo, eleValHomo) / xy[0];
			T_2[2 * i] = T_1[2 * i] + ep * ep * corrector21;

			corrector22 = MVO.dot_product(N11H11[1], eleIdxCell, eleValCell)
					* MVO.dot_product(SeU[0], eleIdxHomo, eleValHomo);
			corrector22 += MVO.dot_product(F11R11[1], eleIdxCell, eleValCell)
					* MVO.dot_product(SePhi[0], eleIdxHomo, eleValHomo);
			corrector22 += MVO.dot_product(PW[1], eleIdxCell, eleValCell);
			corrector22 += MVO.dot_product(NN1HH1[1], eleIdxCell, eleValCell)
					* MVO.dot_product(gradU[0], eleIdxHomo, eleValHomo) / xy[0];
			corrector22 += MVO.dot_product(MM1SS1[1], eleIdxCell, eleValCell) * uVal / (xy[0] * xy[0]);
			corrector22 += MVO.dot_product(FF1RR1[1], eleIdxCell, eleValCell)
					* MVO.dot_product(gradPhi[0], eleIdxHomo, eleValHomo) / xy[0];
			T_2[2 * i + 1] = T_1[2 * i + 1] + ep * ep * corrector22;
		}
	}

	public void assemblePEEigenSOTSSolution(FES1L2PE fsFine, double[][] v0, double[][] v1, double[][] v2,
                                            ScalarFunc rhoFunc, double[] homoD, double[] d1, double[] d2, double ep) {
		// TODO Auto-generated method stub
		int numEigen = v0.length;
		int nv = fsFine.mesh.getNv();
		FES1L21 fs = new FES1L21((Mesh1L2)fsFine.mesh);
		double[] rhoVal = fs.valueFromFunc(rhoFunc, null);
		double[][] v0_uphi = new double[2][nv];
		double[][] v1_uphi = new double[2][nv];
		double[][] v2_uphi = new double[2][nv];
		double[] u0u0phi = new double[nv];
		double[] u0u1phi = new double[nv];
		double[] u0u2_u1u1phi = new double[nv];
		for(int i = 0; i < numEigen; i++) {
			fsFine.extractComponent(v0[i], v0_uphi);
			fsFine.extractComponent(v1[i], v1_uphi);
			fsFine.extractComponent(v2[i], v2_uphi);
			for(int j = 0; j < nv; j++) {
				u0u0phi[j] = rhoVal[j] * v0_uphi[0][j] * v0_uphi[0][j];
				u0u1phi[j] = rhoVal[j] * (v1_uphi[0][j] - v0_uphi[0][j]);
				u0u2_u1u1phi[j] = rhoVal[j] * ((v2_uphi[0][j] - v0_uphi[0][j]) * v0_uphi[0][j] 
						                      + (v1_uphi[0][j] - v0_uphi[0][j]) * (v1_uphi[0][j] - v0_uphi[0][j]));
			}
			double lambda1 = - homoD[i] * fs.integrate(u0u1phi) / fs.integrate(u0u0phi) / ep;
			d1[i] = homoD[i] + ep * lambda1;
			d2[i] = d1[i] -homoD[i] * fs.integrate(u0u2_u1u1phi) / fs.integrate(u0u0phi) 
					+ ep * ep * 2 * lambda1 * lambda1 / homoD[i];
		}
		
	}
	
	public void assembleAxisPEEigenSOTSSolution(FES1L2PE fsFine, double[][] v0, double[][] v1, double[][] v2,
                                                ScalarFunc rhoFunc, double[] homoD, double[] d1, double[] d2, double ep) {
		// TODO Auto-generated method stub
		int numEigen = v0.length;
		int nv = fsFine.mesh.getNv();
		FES1L21 fs = new FES1L21((Mesh1L2)fsFine.mesh);
		double[] rhoVal = fs.valueFromFunc(rhoFunc, null);
		double[][] v0_uphi = new double[2][nv];
		double[][] v1_uphi = new double[2][nv];
		double[][] v2_uphi = new double[2][nv];
		double[] u0u0phi = new double[nv];
		double[] u0u1phi = new double[nv];
		double[] u0u2_u1u1phi = new double[nv];
		double[] x = new double[1];
		for(int i = 0; i < numEigen; i++) {
			fsFine.extractComponent(v0[i], v0_uphi);
			fsFine.extractComponent(v1[i], v1_uphi);
			fsFine.extractComponent(v2[i], v2_uphi);
			for(int j = 0; j < nv; j++) {
				fsFine.mesh.getCoordPerNode(j, x);
				u0u0phi[j] = rhoVal[j] * v0_uphi[0][j] * v0_uphi[0][j] * x[0];
				u0u1phi[j] = rhoVal[j] * (v1_uphi[0][j] - v0_uphi[0][j]) * x[0];
				u0u2_u1u1phi[j] = rhoVal[j] * x[0] * ((v2_uphi[0][j] - v1_uphi[0][j]) * v0_uphi[0][j] 
						                      + (v1_uphi[0][j] - v0_uphi[0][j]) * (v1_uphi[0][j] - v0_uphi[0][j]));
			}
			double lambda1 = - homoD[i] * fs.integrate(u0u1phi) / fs.integrate(u0u0phi) / ep;
			d1[i] = homoD[i] + ep * lambda1;
			d2[i] = d1[i] - homoD[i] * fs.integrate(u0u2_u1u1phi) / fs.integrate(u0u0phi) 
					+ ep * ep * 2 * lambda1 * lambda1 / homoD[i];
		}
		
	}
	
	FES1L2PE fsPE;
}
